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INFLUENCE OF FINITE ELEMENT SOFTWARE ON ENERGY RELEASE RATES 
COMPUTED USING THE VIRTUAL CRACK CLOSURE TECHNIQUE 

Ronald Krueger and Dirk Goetze + 


ABSTRACT 

Strain energy release rates were computed along straight delamination fronts of Double 
Cantilever Beam, End-Notched Flexure and Single Leg Bending specimens using the 
Virtual Crack Closure Technique (VCCT). The results were based on finite element 
analyses using ABAQUS” and ANSYS ' and were calculated from the finite element 
results using the same post-processing routine to assure a consistent procedure. Mixed- 
mode strain energy release rates obtained from post-processing finite element results 
were in good agreement for all element types used and all specimens modeled. 
Compared to previous studies, the models made of solid twenty-node hexahedral 
elements and solid eight-node incompatible mode elements yielded excellent results. For 
both codes, models made of standard brick elements and elements with reduced 
integration did not correctly capture the distribution of the energy release rate across the 
width of the specimens for the models chosen. The results suggested that element types 
with similar formulation yield matching results independent of the finite element 
software used. For comparison, mixed-mode strain energy release rates were also 
calculated within ABAQUS ^/Standard using the VCCT for ABAQUS" add on. For all 
specimens modeled, mixed-mode strain energy release rates obtained from ABAQUS” 
finite element results using post-processing were almost identical to results calculated 
using the VCCT for A BAQUS* add on. 


1. INTRODUCTION 

One of the most common failure modes for composite structures is delamination [1-4]. The 
remote loadings applied to composite components are typically resolved into interlaminar 
tension and shear stresses at discontinuities that create mixed-mode I, II and III delaminations. 
To characterize the onset and growth of these delaminations, the use of fracture mechanics has 
become common practice over the past two decades [4-6], The total strain energy release rate, 
Gt, the mode I component due to interlaminar tension, Gi, the mode II component due to 
interlaminar sliding shear, Gn, and the mode III component, Gm, due to interlaminar scissoring 
shear, as shown in Figure 1, have to be calculated when a methodology based on fracture 
mechanics is used. In order to predict delamination onset or growth for two-dimensional 
problems, these calculated G components are compared to interlaminar fracture toughness 
properties measured over a range from pure mode I loading to pure mode II loading [7-9]. A 
quasi static mixed-mode fracture criterion is determined by plotting the interlaminar fracture 
toughness, G c , versus the mixed-mode ratio, Gu/Gt, determined from data generated using 
Double Cantilever Beam (DCB) (Gn/Gj=0), End Notched Flexure (ENF) (Gn/Gf =\ ), and Mixed 
Mode Bending (MMB) tests of varying ratios, as shown in Figure 2 for T300/914C [10]. A curve 
fit of these data is performed to determine a mathematical relationship between G c and Gi/Gt- 
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[11]. Failure is expected when, for a given mixed-mode ratio Gu IGt, the calculated total energy 
release rate, Gt, exceeds the interlaminar fracture toughness, G c . An interaction criterion 
incorporating the scissoring shear (mode III), however, has not yet been established. The edge- 
cracked torsion test (ECT) is being considered for standardization [12, 13]. 

The virtual crack closure technique (VCCT) is widely [14, 15] used for computing energy 
release rates based on results from continuum (2D) and solid (3D) finite element analyses to 
supply the mode separation required when using the mixed-mode fracture criterion. Although the 
original publication on VCCT dates back a quarter century [14], the virtual crack closure 
technique has only recently been implemented into ABAQUS® [16]. Currently the virtual crack 
closure technique (VCCT) is used mainly by scientists in universities, research institutions and 
government laboratories and is usually implemented in specialized codes or used in post- 
processing routines in conjunction with general-purpose finite element software. 

Lately, an increased interest in using a fracture mechanics based approach to assess the 
damage tolerance of composite structures in the design phase and during certification has also 
renewed the interest in the virtual crack closure technique. Reliability and repeatability of results 
are required to ensure that methodologies based on fracture mechanics will be incorporated into 
the design and certification process. Computed energy release rates as shown in Figure 3 for the 
two-dimensional finite element analysis of a hybrid composite tapered flex-beam have shown 
that energy release rates computed from results obtained from different finite element codes are 
in good agreement [17]. 

The first objective of this study was to compare computed strain energy release rates using 
VCCT based on results obtained from three-dimensional analyses using different finite element 
codes. For comparison, the finite element codes ANSYS ' and ABAQUS® 1 were chosen which 
were run on different platforms. For this investigation, simple specimens like the double 
cantilever beam (DCB) specimen and the End-Notched Flexure (ENF) specimen with 
unidirectional layup and the Single Leg Bending (SLB) specimen with multi- directional layup 
were chosen for full three-dimensional finite element simulations. The specimens were modeled 
with different types of solid brick elements. Mixed-mode strain energy release rates were 
calculated using a user written post-processing routine. Computed strain energy release rates 
across the width of the specimens were compared with results from previous analyses [18, 19]. 

The second objective of this study was to test and study the implementation of the virtual 
crack closure technique in ABAQUS ® [16]. Analyses using models of DCB and ENF specimens 
with unidirectional layup and models of a SLB specimen with multi-directional layup were 
performed. Computed mixed-mode strain energy release rates across the width of the specimens 
were compared with results obtained from a user written post-processing routine and with results 
from previous analyses [18, 19]. 


2. SPECIMEN DESCRIPTION 

For the current numerical investigation simple specimens such as the Double Cantilever 
Beam (DCB), the End-Notched Flexure (ENF) and the Single Leg Bending (SLB) specimens, as 
shown in Figure 4, were chosen. The DCB specimen is used to determine the mode I 
interlaminar fracture toughness, Gic (Gn/Gj=0), the ENF specimen is used for mode II 
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(Gn/G]=l) fracture toughness testing. The SLB specimen was introduced for the determination 
of fracture toughness as a function of mixed-mode I/II ratio [20, 21]. This test may be performed 
in a standard three point bending fixture such as that used for the ENF test. By varying the 
relative thickness of the delaminated regions (ti and t 2 ) various mode mixities may be achieved. 
This type of specimen was chosen to study mode separation. Previously, a number of combined 
experimental and numerical studies of these specimens had been performed where the critical 
strain energy release rates were evaluated under mode I, mode II and mixed-mode conditions 
[18-20, 22, 23], 

In general, mode I, mode II and mixed-mode tests are performed on unidirectionally 
reinforced laminates, which means that delamination growth occurs at a [0/0] interface and crack 
propagation is parallel to the fibers. For the current study, DCB and ENF specimens made of 
T300/1076 graphite/epoxy and unidirectional layup [0] 24 (UD24) were modeled. One SLB- 
specimen modeled was made of C12K/R6376 graphite/epoxy and unidirectional layup [0] 32 
(UD32). Although this unidirectional layup is desired for standard test methods to characterize 
fracture toughness, this kind of delamination growth will rarely occur in real structures. 
Previously, a number of combined experimental and numerical studies on specimens with 
multidirectional layup were performed where the critical strain energy release rates of various 
interfaces were evaluated under mode I, mode II and mixed-mode conditions [20, 22, 23]. 
Therefore, another SFB-specimen made of C12K/R6376 graphite/epoxy with a multidirectional 
layup was also modeled. The stacking sequence [±30/0/-30/0/30/0 4 /30/0/-30/0/-30/30/ t - 
30/30/0/30/0/ -30/0 4 /30/0/30/0/±30] was designated D±30, where the arrow denotes the location 
of the delamination. The material properties are given in Table 1 and the layups modeled are 
summarized in Table 2. 


3. ANALYSIS TOOLS 

3.1 Virtual Crack Closure Technique 

For delaminations in laminated composite materials where the failure criterion is highly 
dependent on the mixed-mode ratio (as shown in Figure 2), the virtual crack closure technique 
(VCCT) has been widely used for computing energy release rates [14, 15]. Results based on 
continuum (2-D) and solid (3-D) finite element analyses provide the mode separation required 
when using the mixed-mode fracture criterion. 

The mode I and mode II components of the strain energy release rate, G„ G u are computed 
as shown in Figure 5(a) for a 2-D four-node element as an example of VCCT. The terms X\ , Z\ 
are the forces at the crack tip at nodal point i and u’ £ , w’ £ («’,*, w’ p ) are the displacements at the 
corresponding nodal points t and t behind the crack tip. For geometric nonlinear analysis where 
large deformations may occur, both forces and displacements obtained in the global coordinate 
system need to be transformed into a local coordinate system (x\ z") which originates at the crack 
tip as shown in Figure 5(a). The local crack tip system defines the tangential ( x' or mode II) and 
normal (?' or mode I) coordinate directions at the crack tip in the deformed configuration. The 
total energy release rate G r is calculated from the individual mode components as G r =G, +G U 
+G m , where G IU =0 for the two-dimensional case shown in Figure 5(a). The extension to 3-D is 
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straight forward as shown in Figure 5(b). Equations and additional details about VCCT for 
higher order solid elements and shell elements can be found in the literature [14, 15]. 


3.2 User Written VCCT Post-Processing Routines 

Although the original publication on VCCT dates back a quarter century [14], the virtual 
crack closure technique has only recently been implemented into ABAQUS ' [16]. Currently, the 
virtual crack closure technique (VCCT) is used mainly by scientists in universities, research 
institutions, government laboratories and research groups in industry who use their own user 
written subroutines. The subroutines may either interface directly with the finite element code 
during its execution, provided this option has been made available for this particular software, or 
operate as an entirely separate post-processing step. 


3.2.1 Post-Processing Routine for ABAQUS® 

The application of the virtual crack closure technique based on results from finite element 
analysis requires access to the element forces at nodes, the nodal point displacements and the 
nodal point coordinates. The flow chart, depicted in Figure 6, shows an independent post- 
processing routine where input data for VCCT is extracted directly from an ABAQUS B binary 
result file (.fil): 

1. Establish interface with finite element software ABAQUS ®. 

2. Provide input interface to read variable problem- specific external data such as 
crack length, Young’s modulus, etc., and control parameters such as element type, 
etc., which is generally not hard coded into the subroutine. 

3. Create output file to store echo print of control parameters, input data, and retrieved 
data as well as intermediate and final results. 

4. Read element forces, or forces at constraints, nodal point coordinates, and nodal 
displacements from binary result file (.fil). 

5. Store the data retrieved in step 4 in internal common blocks and arrays for further 
access during the calculation. 

6. Calculate area virtually closed, A A, using nodal point coordinates. 

7. Calculate local coordinate system x’, y’ and z’ at crack tip for geometrically 
nonlinear analysis and arbitrarily shaped delamination front as discussed in section 
3.1 or related references [14, 15]. 

8 . 

a. Obtain the forces at the crack tip and in front of the crack tip from the 
forces at element nodes by summing the forces at common nodes from 
elements belonging to the upper surface. 
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b. Obtain relative displacements between the surfaces from the nodal point 
displacements. 

c. Transform forces and displacements to local crack tip coordinate system. 

9. Calculate mixed-mode strain energy release rates using the equations given for the 
individual element type and correct for different element lengths and widths if 
necessary. 

10. Repeat steps 6 through 9 for all nodes along the delamination front and write 
results to output fde. 


3.2.2 Post-Processing Routine for ANSYS Mechanics™ 

The basic strategy to implement VCCT in a post-processing routine is the same regardless 
of data retrieved from a binary data file such as the ABAQUS® .fil file or a data file in ASCII 
format such as a PATRAN® 1 neutral file. The original post-processing routine to calculate energy 
release rates was modified to allow more flexibility. It is finite element software independent and 
uses three ASCII files as input which need to be provided by the user. The modified approach 
shown in Figure 7 can be used with any finite element software as long as the input data files 
required can be created. It was tested in conjunction with the finite element software ANSYS® 1 , 
for which a macro ( vcct-dat.mac ) was written to create the data files required. In a first step, the 
model data was written to the file model.dat. The nodal coordinates for elements at the 
delamination front are written to the file topbot.dat, the displacements and the forces are written 
to the file residt.dat. 

In the post-processing routine, the software specific interface to the binary output file in 
step 1 of the procedure became obsolete as shown in Figure 7. Major changes were also made to 
allow the routine to interface with the ASCII data files model.dat, topbot.dat and residt.dat and 
to read the data. Additional changes were also required to allow for the conversion of special 
element designations unique to ANSYS® to a generic designation used in the user written post- 
processing routine. 


3.3 VCCT for ABAQUS® 

Currently, VCCT for ABAQUS ® is an add-on capability to ABAQUS ®/Standard Versions 
6.5 and 6.6 that provides a specific implementation of the virtual crack closure technique within 
ABAQUS®. This new capability enhances ABAQUS® to solve delamination and debonding 
problems in composite materials. The implementation is compatible with all the features in 
ABAQUS 8 such as large-scale nonlinear, models of composite structures including continuum 
shells, composite materials, cohesive elements, buckling, and contact. The bond line in two- 
dimensional analyses or the plane of delamination in three-dimensional analyses is modeled using 
the existing ABAQUS ®/Standard crack propagation capability based on the contact pair 
capability. Additional element definitions are not required, and the underlying finite element mesh 
and model does not have to be modified [16]. 

Beyond simple calculation of the mixed-mode strain energy release rates at the crack tip or 
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along the delamination front, the implementation also offers a crack propagation capability in 
ABAQUS* . It is implied that the energy release rate at the crack tip is calculated at the end of a 
converged increment. Once the energy release rate exceeds the critical strain energy release rate 
(including the user- specified mixed-mode criteria as shown in Figure 2), the node at the crack tip 
is released in the following increment, which allows the crack to advance. To avoid sudden loss of 
stability when the crack tip is advanced, the force at the crack tip before advance is released 
gradually during succeeding increments in such a way that the force is brought to zero no later 
than the time at which the next node along the crack path begins to open. The propagation 
capability was not tested since it exceeded the scope of the current study. 


4 FINITE ELEMENT MODELING AND ANALYSIS 
4.1 Analysis Procedure 

The goal of this investigation was to compare computed strain energy release rates based 
on results obtained from different finite element codes. For comparison, the finite element codes 
ANSYS® and ABAQUS® were chosen. The input files for both codes were generated by a user 
written FORTRAN routine (genmesh. f) as shown in the flow chart of Figure 8. This approach 
was used to avoid user interference with input data and ensure that nodal point coordinates, 
element topology and material data as well as load and boundary conditions of the models were 
identical for both codes. Results from ABAQUS® were saved in the binary .fil results file, results 
from ANSYS” were written to three output files (model.dat, results.dat and 
topbot . dat) using a user written ANSYS” macro as mentioned above. For both codes, mixed- 
mode strain energy release rates were calculated from finite element results using the same user 
written post-processing routine for VCCT to assure a consistent procedure. For comparison, 
mixed-mode strain energy release rates were also calculated within ABAQUS “/Standard using 
the VCCT for ABAQU^ add on mentioned above. 


4.2 Finite Element Selection 

To study the influence of element selection on the computed mixed-mode strain energy 
release rates, several three-dimensional solid element types were used in the past to model the 
specimens [18, 19]. In ABAQUS®, the use of standard solid eight-node brick element C3D8, 
incompatible mode element C3D8I, reduced integration element C3D8R as well as solid twenty- 
node hexahedral elements C3D20 and reduced integration element C3D20R was studied [19]. 
Shear locking is common in first-order, fully integrated elements, such as C3D8, that are subject 
to bending. The numerical formulation of this element gives rise to shear strains that do not 
really exist. Therefore, these elements are too stiff in bending, and many elements over the 
thickness are required to obtain acceptable results. Elements where a lower-order, reduced 
integration is used to form the element stiffness such as the C3D8R and C3D20R elements 
usually provide more accurate results in bending and reduce running time. Incompatible mode 
elements, such as C3D8I are also recommended for bending and contact problems. In these 
elements, internal deformation modes are added to the standard displacement modes of the 
element in order to eliminate the parasitic shear stresses that occur in bending [24], The use of 
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continuum shell elements was also investigated. The continuum shell elements SC8R provided in 
ABAQUS® allow the user to discretize the entire three-dimensional body. The elements look like 
three-dimensional solids from a modeling standpoint and have only displacement degrees of 
freedom. The continuum shell elements can be stacked to provide more refined through- 
thickness response. Due to the kinematic and constitutive formulation of the element and the 
first-order layer-wise composite theory used for layered elements, the continuum shell elements 
appear attractive for modeling [24]. 

Three different eight-node elements were chosen from the ANSYS " element library [25, 
26]. The solid elements SOLID45, SOLID46, and SOLID185 were used for three-dimensional 
modeling of the specimen. Analyses with elements SOLID45 and SOLID46 can be performed 
with extra displacement shapes. The element formulation introduces extra internal degrees of 
freedom (inaccessible to ANSYS® users) to modify the displacement/velocity gradient tensor. 
All internal degrees of freedom are introduced automatically at the element level and are 
condensed. Element SOLID46 additionally allows up to 250 different material layers. Element 
SOLID 185 does not have extra shapes and therefore can exhibit shear locking in bending types 
of applications. In such situations, the use of a uniform reduced integration technique with 
hourglass stiffness for controlling hourglass modes is recommended (SOLID 185R). which is also 
available for element SOLID45 (SOLID45R). The analysis, however, cannot capture the bending 
behavior with a single layer of elements. For example, in the case of a fixed-end cantilever with 
a lateral point load four elements over the thickness are usually recommended. Element 185 also 
has an enhanced strain mode (185E). This mode prevents shear locking in bending-dominated 
problems and volumetric locking in nearly incompressible cases. 

Two twenty-node elements were chosen from the ANSYS® element library. The solid 
elements used in this study are 95 and 186 which are higher order versions of the 45 and 185 
types mentioned above and exhibit quadratic displacement behavior. The elements can tolerate 
irregular shapes without as much loss of accuracy and are well suited to model curved 
boundaries. Element 186 is available as a structural solid and layered solid [25, 26], 


4.3 Finite Element Models 

Typical three-dimensional finite element models of double cantilever beam (DCB), End- 
Notched Flexure (ENF) and Single Leg Bending (SLB) specimens are shown in Figures 9 
through 14. For the DCB and SLB specimen made with unidirectional layup only half of the 
specimen width (B/2) was modeled as shown in Figure 9, and symmetry conditions were applied 
in the xz-plane (y=0). Along the length, all models were divided into different sections with 
different mesh refinement as shown in Figure 9(a). A refined mesh of length d= 5 mm with 20 
elements as used for the UD24 layup is shown in the detail of Figure 9(b), This section length, d, 
had been selected in previous studies [18, 19] and was also used during the current 
investigations. Across the width, all models were divided into a center section of width, i, and a 
refined edge section, j, to capture local edge effects and steep gradients. These sections appear as 
dark areas in the full view of the specimen as shown in Figure 9. The specimens with 
unidirectional layup were modeled with six elements through the specimen thickness ( 2h ) as 
shown in the detail of Figure 9(b). The delamination was modeled as a discrete discontinuity in 
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the center of the specimen, with separate, unconnected nodes (with identical coordinates) on the 
upper and lower surfaces of the delaminated section. For the DCB with unidirectional layup, 
finite element meshes were also developed to model the full width B as shown in Figures 10 and 
1 1. A model with a refined mesh near the edges is shown in Figure 10(a) (/= 1.0 mm, 8 elements) 
and a model with uniform mesh across the width is shown in Figure 10(b). The difference in the 
model shown in Figure 1 1 is a refined zone along the length near the delamination front (<i=10 
mm, 20 elements). 

For the analysis of the ENF specimen, a model with refined edges was used as shown in 
Figure 12. Interpenetration of the cantilever arms was prevented by introducing multi -point 
constraints in the plane of delamination along a line of nodes above the support as discussed in a 
previous study [19]. The use of multi-point constraints appears advantageous as less modeling 
effort is required, and a computationally expensive contact analysis is avoided. 

For the SLB specimen, a length of d= 3.0 mm with 12 elements in the vicinity of the 
delamination front was chosen as shown in Figure 13. This section length, d, had been selected 
in previous studies [18] and remained unchanged during the current investigations. For the 
specimen with D±30 layup, two plies on each side of the delamination were modeled 
individually using one element for each ply as shown in the detail of Figure 13. The adjacent 
four plies were modeled by one element with material properties smeared using the rule of 
mixtures [27, 28], This procedure did not calculate the full A-B-D stiffness matrix contributions 
of the plies [27, 28]; however, it appeared suitable to enforce a reasonable model size. The 
adjacent element extended over the four 0° plies. The six outermost plies were modeled by one 
element with smeared material properties. Additionally, two models with a finer mesh across the 
width were created a shown in Figure 14. A model with a refined mesh near the edges is shown 
in Figure 14(a) (j= 2.0 mm, 8 elements) and a model with uniform mesh across the width is 
shown in Figure 14(b). 


5. COMPUTED MIXED-MODE ENERGY RELEASE RATES 

5.1 Results Obtained from ABAQUS* and ANSYS® Using Post-Processing 

The first objective of this study was to compare computed strain energy release rates using 
VCCT based on results obtained from three-dimensional analyses using different finite element 
codes. For comparison, the finite element codes ANSYS 8 and ABAQUS 8 were chosen. For this 
investigation, double cantilever beam (DCB) specimens and the Single Leg Bending (SLB) 
specimens with unidirectional and multi-directional layup were chosen for full three-dimensional 
finite element simulations [26]. Mixed-mode strain energy release rates were calculated using a 
user written post-processing routine as discussed above. Computed strain energy release rates 
across the width of the specimens were compared with results from previous analyses [18, 19]. 

5.1.1 Mode I Energy Release Rates Distribution Across the Width of a DCB Specimen 

The computed mode I strain energy release rate values across the width of a DCB 
specimen with unidirectional layup were normalized with respect to the value from beam theory, 



I, beam 


(«) = 


l2-a 2 -P 2 
B 2 • h 2 ■ E l 


where a denotes the delamination length, P the external load, B the specimen width, h the 
thickness of the cantilever arms as shown in Figure 4 and Ei the modulus of elasticity [8]. The 
normalized values were plotted versus the normalized width, y/B, of the specimen as shown in 
Figure 15 for results obtained from ABAQUS ' . Strain energy release rates obtained from the 
model shown in Figure 9 using C3D8I, C3D20 and C3D20R elements are in excellent agreement 
with the values from previous analyses, which were added for comparison [18], The mode I 
strain energy release rate is fairly constant in the center part of the specimen progressively 
decreasing towards the edges as explained earlier which will cause the initial straight front to 
grow into a curved front as explained in detail in the literature [29-32], As expected, the mode II 
and mode III strain energy release rates were computed to be nearly zero and hence are not 
shown. The model made of C3D8 elements yielded results, which did not correctly capture the 
decrease of the mode I strain energy release rate towards the edge. The mode I strain energy 
release rate computed from the model made of C3D8R elements was noticeably higher across the 
entire width compared to the other results. 

The results of a brief mesh refinement study for the standard eight-node brick element 
(C3D8) and eight node brick elements with reduced integration and hourglass control (C3D8R) 
are shown in Figures 16 to 18 for a uniform mesh across the width of the specimen. Two 
elements over the thickness of one arm were used to compute the distributions shown in 
Figure 16. Three elements over the thickness of one arm yielded the results in Figure 17, and 
four elements over the thickness of one arm yielded the results shown in Figure 18. With 
increasing number of elements, the strain energy release rates distributions approach the 
reference results. Hence, a minimum number of four elements across the thickness is 
recommended when using these elements. Results obtained from models using eight node 
incompatible mode elements were included in the plots of Figures 16 to 18. The distributions 
suggest that as few as two elements across the thickness are sufficient to obtain a good 
agreement with reference results. 

Mode I energy release rates obtained from the model shown in Figure 9 using the finite 
element software ANSYS® are shown in Figure 19 for all element types used in this study. Strain 
energy release rates obtained using SOLID 186, SOLID95, SOLID45, SOLID46 and 
SOLID 185E type elements were in excellent agreement with the values from previous analyses, 
which were added for comparison [18]. The model made of SOLID 185 elements yielded results, 
which did not correctly capture the decrease of the mode I strain energy release rate towards the 
edge. The mode I strain energy release rate computed from the model made of SOLID45R and 
SOLID 185R elements was noticeably higher across the entire width compared to the other 
results. 

Strain energy release rates obtained from analyses using ABAQUS 8 and ANSYS' 8 ' are 
plotted in Figures 20 to 23 for corresponding element types and compared to results obtained 
from previous analyses [18]. For both finite element codes, mode I energy release rates obtained 
from the models using higher order twenty-node elements are in excellent agreement with the 
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values from previous analyses as shown in Figure 20. The models made of standard eight node 
brick elements yielded results that did not correctly capture the decrease of the mode I strain 
energy release rate towards the edge as shown in Figure 21. The mode I strain energy release rate 
computed from models made of eight node brick elements with reduced integration were 
noticeably higher across the entire width compared to the other results obtained from previous 
analyses [18] as shown in Figure 22. For both finite element codes, mode I energy release rates 
obtained from the models using solid eight node incompatible mode elements are in excellent 
agreement with the values from previous analyses as shown in Figure 23. The results obtained 
from the model shown in Figure 9 suggest that element types with similar formulation yield 
matching results independent of the software used. Further, mode I energy release rates obtained 
from the models using higher order twenty-node elements and eight node incompatible mode 
elements are in excellent agreement with the values from previous analyses. 

5,1.2 Mixed-Mode Energy Release Rate Distributions Across the Width of a Unidirectional 
SLB Specimen 

For this investigation, the symmetry of the SLB specimen was taken into account, and 
only one half of the specimen width B/2 was modeled as shown in Figure 9. The finite element 
model is basically identical to the one used for the DCB specimens discussed in the previous 
section, except for the boundary conditions. The influence of element selection on the computed 
mixed-mode strain energy release rates was studied using ABAQUS® C3D8, C3D8I, C3D8R, 
C3D20 and C3D20R type elements as well as ANSYS® SOLID45 and 45R, SOLID 185, 185R 
and 185E and SOLI46, SOLID95 and SOLID 186 type elements [26]. 

In Figures 24 to 26, the computed mode I, II and mode III strain energy release rates from 
ABAQUS® analyses are plotted versus the normalized width, y/B, of the specimen. Computed 
strain energy release rates obtained from models using C3D8I, C3D20 and C3D20R elements are 
in excellent agreement with the values from previous analyses [18]. As shown in Figure 24, the 
mode I strain energy release rate is fairly constant in the center part of the specimen and 
progressively decreases towards the edges as previously discussed for the DCB specimen. The 
model made of C3D8 elements yields results that do not correctly capture the decrease of the 
mode I strain energy release rate towards the edge. The values computed from the model made 
of C3D8R elements appear excessively high. The mode II strain energy release rate as shown in 
Figure 25 is fairly constant across almost the entire width of the specimen, peaking in the 
immediate vicinity of the edges. As shown in Figure 26, the mode III contribution is zero in the 
center of the specimen peaking to about only 8% of Gu at the edges. The model made of C3D8 
elements yields results that do not correctly capture the mode II and III distribution. The mode II 
values computed from the model made of C3D8R elements appear excessively high across the 
entire width of the specimen. These observations support the results obtained from the study of 
the UD24 DCB specimen discussed above. Compared to the model made of C3D20 elements, 
the models made of C3D8I elements yield nearly the same results and provide a reduced model 
size. 

Computed mode I, II and mode III strain energy release rates from ANSYS® analyses are 
plotted in Figures 27 to 29 versus the normalized width, y/B, of the specimen. Computed strain 
energy release rates obtained from models using SOLID 186, SOLID95, elements are in excellent 
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agreement with the values from a previous analysis models [18] as shown in Figure 27. The 
model made of SOLID 185 elements yields results that do not correctly capture the decrease of 
the mode I strain energy release rate towards the edge. The values computed from the model 
made of SOLID45R and SOLID 185R elements appear excessively high. The mode II strain 
energy release rate, as shown in Figure 28, is fairly constant across almost the entire width of the 
specimen, peaking in the immediate vicinity of the edges. As shown in Figure 29, the mode III 
contribution is zero in the center of the specimen peaking to about only 8% of Gn at the edges. 
The model made of SOLID 185 elements yields results that do not correctly capture the mode II 
and III distribution. The mode II values computed from the model made of SOLID45R and 
SOLID 185R elements appear excessively high across the entire width of the specimen. These 
observations support the results obtained from the study of the UD24 DCB specimen discussed 
above. Compared to the models made of SOLID 186, SOLID95 type elements, the models made 
of SOLID45, SOLID46 and SOLID 185E type elements yield nearly the same results that are also 
in excellent agreement with reference results [18]. 

Mixed-mode strain energy release rates obtained from analyses using ABAQUS® and 
ANSYS® are plotted in Figures 30 to 41 for corresponding element types and compared to 
results obtained from previous analyses [18]. For both finite element codes, mode I, II and III 
energy release rates obtained from the models using higher order twenty-node elements are in 
excellent agreement with the values from previous analyses as shown in Figures 30 to 32. The 
models made of standard eight node brick elements yielded results that did not correctly capture 
the distribution of the energy release rates across the width as shown in Figure 33 to 35 for the 
individual fracture modes. The strain energy release rates computed from models made of eight 
node brick elements with reduced integration were noticeably higher across the entire width 
compared to the other results obtained from previous analyses [18] as shown in Figure 36 to 38 
for the individual fracture modes. For both finite element codes, mixed-mode energy release 
rates obtained from the models using solid eight node incompatible mode elements are in 
excellent agreement with the values from previous analyses as shown for the individual models 
in Figures 39 to 41. The results suggest that element types with similar formulation yield 
matching results independent of the finite element software used. Further energy release rates 
obtained from the models using higher order twenty-node elements and eight node incompatible 
mode elements are in excellent agreement with the values from previous analyses. Therefore, 
solid eight node incompatible mode type elements (ABAQUS® C3D8I and ANSYS® SOLID46) 
were chosen for the following study. 

5.1.3 Mixed-Mode Energy Release Rate Distributions Across the Width of a 
Multidirectional SLB Specimen 

In this study, SLB specimens with a multidirectional layup (D±30) were modeled entirely 
with solid elements as shown in Figure 13. The model was made of eight node ABAQUS 8 
C3D8I or ANSYS 8 SOLID46 elements. It was shown above that compared to the models made 
of twenty-node elements, the models made of incompatible mode elements yielded nearly the 
same results and provide a reduced model size. Therefore, these element types were chosen for 
this and the following studies. Section widths and mesh sizes used in the current investigation 
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were taken from a previous study where the effect of mesh size across the width on computed 
mixed-mode energy release had been investigated [18, 19]. 

The computed mode I, II and mode III strain energy release rates, as shown in Figures 42 
to 44, were compared with values from previous analyses [18]. As shown in Figure 42, the zone 
with a constant Gi distribution in the center becomes smaller compared to the UD32 case and the 
decrease towards the edges is more pronounced. The decrease is caused by increased anticlastic 
bending due to the lower values of bending rigidities in the individual arms. Results obtained 
from ABAQUS® and ANSYS " are in excellent agreement and also agree well with the reference 
results. The mode II and mode III distributions are shown in Figures 43 and 44. The mode II 
strain energy release rate is fairly constant across almost the entire width of the specimen and 
peaks near the edges accompanied by local mode III contribution. Compared to the UD32 layup 
these peaks become more visible for specimens with the D±30 layup caused by increased 
anticlastic bending. Results obtained from ABAQUS® follow closely the distributions from the 
reference results. Results obtained from ANSYS® are in good agreement in the center of the 
specimen and also peak towards the edges, however values for y/B< 0 are below and for y/B>0 
above values from ABAQUS 8 and published the reference values [18]. 


5.2 Results Obtained from ABAQUS® Using Post-Processing and VCCT for ABAQUS^ 

The second objective of this study was to test and study the implementation of the virtual 
crack closure technique in ABAQUS® [16]. Analyses using models of DCB and ENF specimens 
with unidirectional layup and models of a SLB specimen with multi- directional layup were 
performed. Computed mixed-mode strain energy release rates across the width of the specimens 
were compared with results obtained from a user written post-processing routine and with results 
from previous analyses [18, 19]. 

5.2.1 Mode I Energy Release Rates Distribution Across the Width of a DCB Specimen 

Strain energy release rates were computed from finite element results obtained from the 
model of a full specimen as shown in Figure 10(a) using C3D8I elements. The incompatible 
mode elements (C3D8I) were chosen since models consistently yielded the best results in the 
studies described above and were least sensitive to mesh refinement. The computed mode I strain 
energy release rate values across the width of a DCB specimen with unidirectional layup were 
normalized with respect to the value from beam theory as above and plotted versus the 
normalized width, y/B, of the specimen as shown in Figure 45. Plotted are results obtained for 
different load increments from VCCT for ABAQUS ® and from the user written post-processing 
routine as shown in the flow chart of Figure 8. Strain energy release rates obtained from both 
methods are practically identical across the entire width and for all load increments. 

Additional analyses were performed using models shown in Figures 10 and 11 which 
allowed a discretization study of the effect of different mesh densities as well as mesh transition 
from coarse in the center to fine at edge. The results plotted in Figure 46 indicate that 
irrespective of mesh density, energy release rates obtained from VCCT for ABAQUS) are almost 
identical to results obtained through post-processing. 


12 



Analyses were repeated with models where the new continuum shell element SC8R was 
used. The computed mode I strain energy release rates are shown in Figure 47 for a uniform 
mesh across the width and a mesh where the edge region had been refined as shown in Figure 10. 
For both models, energy release rates obtained from VCCT for ABAQUS" are almost identical to 
results obtained through post-processing. The results obtained from the model with a uniform 
mesh across the width are in good agreement with previous results, which were included for 
comparison [18]. The model made of SC8R elements however yields different results in the 
refined region near the edge of the specimen when a non-uniform mesh is used. The cause is 
related to the element formulation and the associated problem has been fixed by ABAQUS®. 
After solving the problem, ABAQUS 8 provided a computed mode I strain energy release rate 
distribution obtained from the non-uniform mesh, which is in good agreement with previous 
results. The problem has been fixed for ABAQUS “/Standard 6.6.2 and higher. 

5.2.2 Mode II Energy Release Rates Distribution Across the Width of an ENF Specimen 

For this investigation, the entire ENF specimen was modeled across the width as shown in 
Figure 12. Incompatible mode elements (C3D8I) were chosen since models consistently yielded 
the best results in the studies described above and were least sensitive to mesh refinement. In 
Figure 48, the computed mode II and mode III strain energy release rates are normalized with the 
reference value, Gn.beam, from classical beam theory (not accounting for transverse shear) 

9 • a 2 • P 2 
16- B 2 ■ h 3 • E x ' 

Here, a denotes the delamination length, P the external load, B the specimen width, h the 
thickness of the cantilever arms as shown in Figure 4 and Ej the modulus of elasticity. The 
results are plotted versus the normalized width, y/B, of the specimen. The mode II strain energy 
release rate is fairly constant across almost the entire width of the specimen, peaking in the 
immediate vicinity of the edges. The mode III contribution is zero in the center of the specimen 
peaking to about only 5% of Gn.beam at the edges. The computed G) values are nearly zero and 
therefore are not shown. The results confirm that energy release rates obtained from VCCT for 
ABAQUS 8 are almost identical to results obtained through post-processing. 

5.2.3 Mixed-Mode Energy Release Rate Distributions Across the Width of a 
Multidirectional SLB Specimen 

For the investigation of the SLB specimen with multidirectional layup (D±30), the finite 
element models shown in Figures 13 and 14 were used. The model was made of eight node 
ABAQUS® incompatible mode C3D8I elements. The computed mode I, II and III strain energy 
release rate distributions across the width of the specimen are shown in Figures 49 to 51. The 
zone with a constant G, distribution in the center is smaller compared to the SLB specimen with 
unidirectional layup discussed earlier, and the decrease towards the edges is more pronounced as 
shown in Figure 49. The decrease is caused by increased anticlastic bending due to the lower 
values of bending rigidities in the individual arms. The mode II distribution, as shown in 
Figures 50, is fairly constant across almost the entire width of the specimen and peaks near the 


II, beam 


(«) 
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edges. The mode III contribution is zero in the center of the specimen peaking at the edges as 
shown in Figure 51. Compared to the unidirectional layup (UD32), these peaks become more 
visible for specimens with the multidirectional layup (D±30) caused by increased anticlastic 
bending. The results plotted in Figure 49 to 51 indicate that irrespective of mesh size mixed- 
mode energy release rates obtained from VCCT for ABAQUS” are almost identical to results 
obtained through post-processing. 


6. SUMMARY AND CONCLUDING REMARKS 

Simple fracture toughness specimens such as the Double Cantilever Beam (DCB) 
specimen and the End-Notched Flexure (ENF) specimen with unidirectional layup and the Single 
Leg Bending (SLB) specimen with unidirectional and multi-directional layup were modeled with 
full three-dimensional solid finite elements and analyzed using the finite element codes ANSYS® 
and ABAQUS ®. Model data such as nodal point coordinates, element topology and material data 
as well as load and boundary conditions were generated using a user written routine to guarantee 
that the models were identical for both codes and to avoid interference with input data. Different 
element types were used to model the specimens. Delaminations were modeled as discrete 
discontinuities. Strain energy release rates were calculated along a straight front using the virtual 
crack closure technique. For finite element results based on ABAQUS® and ANSYS®, mixed- 
mode strain energy release rates were calculated from finite element results using the same post- 
processing routine to assure a consistent procedure. For comparison, mixed-mode strain energy 
release rates were also calculated within ABAQUS “/Standard using the VCCT for ABAQU 5® 
add on. 

Mixed-mode strain energy release rates obtained from post processing ANSYS® and 
ABAQUS®’ finite element results were in good agreement with results from previous analyses for 
all specimens modeled. Compared to previous studies, the models made of solid twenty-node 
hexahedral elements (ABAQUS 1 * types C3D20 and C3D20R with reduced integration, ANSYS® 
types SOLID95 and SOLID 186) and solid eight-node incompatible mode elements (ABAQUS® 
type C3D8I, ANSYS® types SOLID45, SOLID46 and SOLID 185E) yielded excellent results. 
The models made of standard brick elements (ABAQUS® type C3D8, ANSYS® type SOLID 185) 
yielded results, which were in agreement but compared to reference results did not correctly 
capture the distribution of the energy release rate across the width of the specimen for the model 
chosen. The values computed from the model made of elements with reduced integration 
(ABAQUS® type C3D8R, ANSYS® types SOLID45R and SOLID 185R) were also in good 
agreement however appeared consistently higher than the reference results. The results suggest 
that element types with similar formulation yield matching results independent of the finite 
element software used. 

For all specimens modeled, mixed-mode strain energy release rates obtained from 
ABAQUS® finite element results using post processing were almost identical to results 
calculated within ABAQUS ®7Standard using the VCCT for ABAQU 5® add on. For identical 
models, almost identical results were computed which indicates that the post-processing routine 
is redundant and may be entirely replaced by the VCCT for ABAQUS* routine. The results also 
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imply that results computed by the VCCT for ABAQUS ® routine can readily be compared to 
existing results computed previously using a post-processing routine. 


ACKNOWLEDGEMENTS 

The authors gratefully acknowledge Mike Sasdelli and Cheryl Liu of ABAQUS 8 East 
and Satish Gandhi of ABAQUS 8 Headquarters for seriously investigating the problem found 
with the continuum shell element and providing a solution for upcoming versions of ABAQUS 8 
Standard. 


REFERENCES 


[1] A.C. Garg, Delamination - A Damage Mode in Composite Structures, Eng. Fracture 
Mech., vol. 29, pp. 557-584, 1988. 

[2] V.V. Bolotin, Delaminations in Composite Structures: Its origin. Buckling, Growth and 
Stability, Composites Part B: Engineering, vol. 27B, pp. 129-145, 1996. 

[3] T.E. Tay, Characterization and Analysis of Delamination Fracture in Composites - An 
Overview of Developments from 1990 to 2001, Applied Mechanics Reviews, vol. 56, pp. 
1-32, 2003. 

[4] T.K. O'Brien, Characterization of Delamination Onset and Growth in a Composite 
Laminate, in Damage in Composite Materials, ASTM STP 775,: American Society for 
Testing and Materials, 1982, pp. 140-167. 

[5] T.K. O'Brien, Fracture Mechanics of Composite Delamination, in ASM Handbook, 
Volume 21, Composites: ASM International, 2001, pp. 241-245. 

[6] R.H. Martin, Incorporating Interlaminar Fracture Mechanics into Design, in 
International Conference on Designing Cost-Effective Composites: IMechE Conference 
Transactions, London, U.K., 1998, pp. 83-92. 

[7] ASTM D 6671-01, Standard Test Method for Mixed Mode 1-Mode 11 Interlaminar 
Fracture Toughness of Unidirectional Fiber Reinforced Polymer Matrix Composites, in 
Annual Book of ASTM Standards, vol. 15.03: American Society for Testing and 
Materials, 2000. 

[8] ASTM D 5528-94a, Standard Test Method for Mode I Interlaminar Fracture Toughness 
of Unidirectional Fiber-Reinforced Polymer Matrix Composites, in Annual Book of 
ASTM Standards, vol. 15.03: American Society for Testing and Materials, 2000. 

[9] T.K. O'Brien, Composite Interlaminar Shear Fracture Toughness, G„ c : Shear 
Measurement or Sheer Myth ?, in Composite Materials: Fatigue and Fracture, Seventh 
Volume, ASTM STP 1330, 1998, pp. 3-18. 

[10] M. Konig, R. Kruger, K. Kussmaul, M. v. Alberti, and M. Gadke, Characterizing Static 
and Fatigue Interlaminar Fracture Behaviour of a First Generation Graphite/ Epoxy 
Composite, in Composite Materials: Testing and Design - (13th Vol.), ASTM STP 1242, 
S .J. Hooper, Ed.: American Society for Testing and Materials, 1997, pp. 60-81. 

[11] M.L. Benzeggagh and M. Kenane, Measurement of Mixed-Mode Delamination Fracture 
Toughness of Unidirectional Glass/Epoxy Composites with Mixed-Mode Bending 
Apparatus, Composites Science and Technology, vol. 56, pp. 439-449, 1996. 


15 



[12] S.M. Lee, An Edge Crack Torsion Method for Mode III Delamination Fracture Testing , 
J. of Composite Technology and Research., pp. 193-201, 1993. 

[13] J.G. Ratcliffe, Characterization of the Edge Crack Torsion (ECT) Test for Mode III 
Fracture Toughness Measurement of Laminated Composites, NASA/TM-2004-2 13269 
September 2004 2004. 

[14] E.F. Rybicki and M.F. Kanninen, A Finite Element Calculation of Stress Intensity 
Factors by a Modified Crack Closure Integral, Eng. Fracture Mech., vol. 9, pp. 931-938, 
1977. 

[15] R. Krueger, Virtual Crack Closure Technique: History, Approach and Applications, 
Applied Mechanics Reviews, vol. 57, pp. 109-143, 2004. 

[16] VCCT for ABAQUS - User's Manual, ABAQUS 2005. 

[17] G.B. Murri, J.R. Schaff, and A.F. Dobyns, Fatigue Life Analysis of Tapered Hybrid 
Composite Flexbeams, presented at ICCM-14, San Diego, 2003. 

[18] R. Kruger, Three Dimensional Finite Element Analysis of Multidirectional Composite 
DCB, SLB and ENF Specimens, Institute for Statics and Dynamics of Aerospace 
Structures, University of Stuttgart ISD-Report No. 94/2, 1994. 

[19] R. Krueger and T.K. O'Brien, A Shell/ 3D Modeling Technique for the Analysis of 
Delaminated Composite Laminates, NAS A/TM-2000-2 10287, ARF-TR-2207, June 
2000 . 

[20] B.D. Davidson, R. Kruger, and M. Konig, Three Dimensional Analysis of Center 
Delaminated Unidirectional and Multidirectional Single Leg Bending Specimens, 
Composites Science and Technology, vol. 54, pp. 385-394, 1995. 

[21] A. Pieracci, B.D. Davidson, and V. Sundararaman, Nonlinear Analyses of Homogeneous, 
Symmetrically Delaminated Single Leg Bending Specimens, Journal Composite Tech. 
Res., vol. 20,’ pp. 170-178, 1998. 

[22] B.D. Davidson, R. Kruger, and M. Konig, Effect of Stacking Sequence on Energy Release 
Rate Distributions in Multidirectional DCB and ENF Specimens, Eng. Fracture Mech., 
vol. 55, pp. 557-569, 1996. 

[23] B.D. Davidson, R. Kruger, and M. Konig, Three Dimensional Analysis and Resulting 
Design Recommendations for Unidirectional and Multidirectional End-Notched Flexure 
Tests, J. Composite Materials, vol. 29, pp. 2108-2133, 1995. 

[24] ABAQUS 6.5, Analysis User's Manual, Volume IV: Elements, ABAQUS 2005. 

[25] ANSYS Elements Reference, ANSYS Inc., Canonsburg, PA 15317 2005. 

[26] D. Goetze, Bruchanalyse von Faserverbunclwerkstoff-Proben mit Finiten Elementen unci 
der Virtuellen R ifischli efiungsmethode ( Computational Fracture Analysis of Composite 
Specimens Based on Three-Dimensional Finite Element Analysis and the Virtual Crack 
Closure Technique), Institut fuer Statik und Dynamik der Fuft- und 
Raumfahrtkonsturktionen, Universitaet Stuttgart, Germany, Diplomarbeit 2003. 

[27] S.W. Tsai, Theory of Composite Design: Think Composites, 1992. 

[28] S.W. Tsai and H.T. Hahn, Introduction to Composite Materials: Technomic Publishing 
Co., Inc., 1980. 

[29] I.S. Raju, K.N. Shivakumar, and J.H. Crews Jr., Three-Dimensional Elastic Analysis of a 
Composite Double Cantilever Beam Specimen, AIAA J., vol. 26, pp. 1493-1498, 1988. 

[30] R. Kruger, M. Konig, and T. Schneider, Computation of Local Energy Release Rates 
Along Straight and Curved Delamination Fronts of Unidirectionally Laminated DCB- 
and ENF - Specimens, in Proceedings of the 34th AIAA/ASME/ASCE/AHS/ASC SSDM 


16 



Conference, La Jolla, CA: American Institute of Aeronautics and Astronautics, 
Washington, 1993, pp. 1332-1342. 

[31] B.D. Davidson and R.A. Schapery, Effect of Finite Width on Deflection and Energy 
Release Rate of an Orthotropic Double Cantilever Specimen, J. Composite Mat., vol. 
22, pp. 640-656, 1988. 

[32] B.D. Davidson, An Analytical Investigation of Delamination Front Curvature in Double 
Cantilever Beam Specimens, J. Composite Mat., vol. 24, pp. 1124-1137, 1990. 


17 



Table 1. 


Material Properties. 


T300/1076 Unidirectional Graphite/Epoxy Prepreg 

E\\ = 139.4 GPa 

£22 = 10.16 GPa 

£33 = 10.16 GPa 

vi 2 = 0.30 

vi 3 = 0.30 

V23 = 0.436 

G 12 = 4.6 GPa 

Gn =4.6 GPa 

G 23 = 3.54 GPa 

C12K/R6376 Unidirectional Graphite/Epoxy Prepreg 

E\\ = 146.9 GPa 

£ 2 2 = 10.6 GPa 

£33 = 10.6 GPa 

vi 2 = 0.33 

vi 3 = 0.33 

V23 = 0.33 

Gi 2 = 5.45 GPa 

G 13 = 5.45 GPa 

G 23 = 3.99 GPa 


The material properties are given with reference to the ply coordinate axes where index 1 1 denotes 
the ply principal axis that coincides with the direction of maximum in-plane Young’s modulus 
(fiber direction). Index 22 denotes the direction transverse to the fiber in the plane of the lamina 
and index 33 the direction perpendicular to the plane of the lamina. 

Table 2. 

Stacking Sequence. 

Layup-ID Stacking Sequence Material 

UD24 [0 ] 24 T300/1076 

UD32 [0 ] 32 C12K/R6376 

D±30 [±30/0/-30/0/30/0 4 /30/0/-30/0/-30/30/ r 30/30/30/0/30/0/-30/0 4 /-30/0/30/0/±30] C12K/R6376 
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Interlaminar sliding shear 
Mode II 

Figure 1: F racture M odes. 



Mode II 



Figure 2. Mixed-mode fracture criterion for T300/914C 
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constant axial load P — 
cyclic displacement v 



(a) Flexbeam specimen with loading and micrograph of observed damage 




delamination length along interface bl, mm 
(b) Computed energy release rate Gji along delaminated interface 
Figure 3. Investigation of hybrid flex-beam specimen [17]. 
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layup UD24 

B 25.0 mm 

2h 3.0 mm 

2L 150.0 mm 

a 1 1 1 .5 mm 

P 12.66 N 


(a) Double Cantilever Beam Specimen (DCB) 



layup UD24 

B 25.0 mm 

2h 3.0 mm 

2L 150.0 mm 

a 30.0 mm 

P 503.0 N 


(b) End Notched Flexure Specimen (ENF) 



layup D±30 

25.4 mm 

2.03 mm 

2.03 mm 
177.8 mm 

34.3 mm 
100.0 N 


Figure 4. Specimen configurations. 
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program grate 

start program 



Figure 6. Flow chart of routine to calculate strain energy release rates from ABAQUS® results 
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program grate 

start program 


vcct-dat .mac 

ANSYS® macro to write 
topology, forces and 



PARAMETER (NRU=1) 


1 © 

read (nunit , *) word, bl , . . . 

read data from user input file 



1 (3) 

call gstart (kunit , output ) 

1 VrV 

open file for output 



call datin (head, 1 , nnum, . . . ) 

read model topology from model, dat 
read elements above the delamination in 
group ELCRTOP, elements below the 
delamination in group ELCRBOT and 
nodes at the delamination front in group 
CRFRONT from topbot. dat 
read displacement at forces at element 
nodes from result.dat 
store data in arrays and common blocks for 
later use 


call eltyp (eltype, elnum, . . .) 

convert ANSYS® element type to internally 
used convention 




call datout (...) 
call orgtopo (nopel, . . . ) 
call varea (iunit , ntype, . 
call local (iunit , ntype, . 
call trans (iunit , ntype) 
call grate (ntype, iunit , . 
calculate area and local coordinate system 
transform forces and displacements, 
compute mixed mode energy release rate 


call gend( iunit) 
end 

close output file and terminate program 


user specified 
input file 
grate.data 
contains parameters 
to choose options 
and input/outputy 
file names 



Figure 7. Flow chart of modified routine to calculate strain energy release rates. 


24 


< user specified 
input file 
genmesh.data 
contains loads, 
dimensions and 
material data 


genmesh.f 

read input data 
calculate nodal point 
coordinates, element topology, 
position of loads and 
boundaries 
write input file 



ABAQUS® Standard 
read input data 

perform finite element analysis 
calculate Gi, Gn, Gm 
write output 


ANSYS Mechanical™ 
read input data 
perform finite element analysis 
write output 



ABAQUS® Viewer 
read output database 
plot results 
write Gi, Gn, Gm 


extract.f 

read data from result file(s) 
calculate Gi, Gn, Gm 
write energy release rates to 
output file 



user specified 
input file 
extract.data 
contains loads 
case, increment 
and 

file names 


Gi, Gn, Gm 


/ Gi, Gn, Gm \ 

I mixed mode ratio ] 

V grate.rpt J 

V failure index J 


x^g rate.txt 


Figure 8. Flow chart of entire analysis procedure. 
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(b) Deformed model ofDCB specimen with uniform mesh 
Figure 10. Full three-dimensional finite element model of a DCB specimen 
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Figure 12. Full three-dimensional finite element model of an ENF specimen 
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— — — — - outline of undeformed configuration 
w c : deflection in the center of the specimen 



Figure 13. Finite element model ofSLB specimen with D+30 layup 









(a) Deformed model ofSLB specimen with refined mesh at the edge 

^ B = i 



2 3 


(b) Deformed model ofSLB specimen with uniform mesh 
Figure 14. Full three-dimensional finite element model of a SLB specimen 
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Energy 
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® J 

Figure 15. Influence of ABAQUS element selection on computed strain energy release rate 
distribution across the width of a DCB specimen with UD24 layup. 
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Figure 16. Influence of element selection on computed strain energy release rate 
distribution across the width of a DCB specimen with UD24 layup. 
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Figure 17. Influence of element selection on computed strain energy release rate 
distribution across the width of a DCB specimen with UD24 layup. 
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Figure 18. Influence of element selection on computed strain energy release rate 
distribution across the width of a DCB specimen with UD24 layup. 
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Figure 19. Influence of AN SYS element selection on computed strain energy release rate 
distribution across the width of a DCB specimen with UD24 layup. 
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Figure 20. Computed strain energy release rate distribution across the width of a 
DCB specimen with UD24 layup modeled with twenty-nodecl solid elements. 
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Figure 21. Computed strain energy release rate distribution across the width of a 
DCB specimen with UD24 layup modeled with eight-noded brick elements. 
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Figure 22. Computed strain energy release rate distribution across the width of a DCB 
specimen with UD24 layup modeled with reduced integration eight-noded brick elements. 
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Figure 23. Computed strain energy release rate distribution across the width of a DCB 
specimen with UD24 layup modeled with incompatible mode eight-nodecl brick elements. 
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Figure 24. Influence of ABAQUS element selection on computed mode I strain energy release rate 
distribution across the width of a SLB specimen with UD32 layup. 
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Figure 25. Influence of ABAQUS® element selection on computed mode II strain energy release 
rate distribution across the width of a SLB specimen with UD32 layup. 
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Figure 26. Influence of ABAQUS® element selection on computed mode 111 strain energy release 
rate distribution across the width of a SLB specimen with UD32 layup. 
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Figure 28. Influence ofANSYS® element selection on computed mode II strain energy release rate 
distribution across the width of a SLB specimen with UD32 layup. 
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Figure 29. Influence of AN SYS' element selection on computed mode 111 strain energy release rate 
distribution across the width of a SLB specimen with UD32 layup. 
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Figure 30. Computed mode I strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with twenty-node solid elements. 
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Figure 31. Computed mode II strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with twenty-node solid elements. 
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Figure 32. Computed mode III strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with twenty-node solid elements. 
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Figure 33. Computed mode I strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with solid eight-node brick elements. 
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Figure 34. Computed mode II strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with solid eight-node brick elements. 
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Figure 35. Computed mode III strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with solid eight-node brick elements. 
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Figure 36. Computed mode I strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with reduced integration eight-node brick elements. 
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Figure 37. Computed mode II strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with reduced integration eight-node brick elements. 
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Figure 38. Computed mode III strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with reduced integration eight-node brick elements. 
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Figure 39. Computed mode I strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with incompatible mode eight-node brick elements. 
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Figure 40. Computed mode II strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with incompatible mode eight-node brick elements. 
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Figure 41. Computed mode I strain energy release rate distribution across the width of a 
SLB specimen with UD32 layup modeled with incompatible mode eight-node brick elements. 
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Figure 42. Influence ofFE software on computed mode I strain energy 

release rate distribution across the width of SLB specimen with D±30 layup. 
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Figure 43. Influence ofFE software on computed mode II strain energy 
release rate distribution across the width ofSLB specimen with D±30 layup. 
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Figure 44. Influence ofFE software on computed mode III strain energy 
release rate distribution across the width ofSLB specimen with D±30 layup. 
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Figure 45. Computed strain energy release rate distribution 
across the width of a DCB specimen for different load increments. 
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Figure 46. Computed strain energy release rate distribution 
across the width of a DCB specimen for different models. 
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Figure 47. Computed strain energy release rate distribution 
across the width of a DCB specimen for continuum shell elements. 
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Figure 48. Computed mixed mode strain energy release rate 
distributions across the width of an ENF specimen. 
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Figure 49. Computed mode I strain energy release rate 
distribution across the width of a SLB specimen for different models. 


Mode II 
Energy 
Release 
Rate, 
N/mm 


0.030 




O 

uniform mesh 


0.025 


X 

VFA 


5 

5 

□ 

refined mesh at the edges 




+ 

VFA 




o 

refined uniform mesh 


0.020 


+ 

VFA 

1 


0.015 - 


0.010 


0.005 - 


0.000 


« * # & ft ® ® • 


-0.4 


- 0.2 


0 

y/B 


0.2 


0.4 


Figure 50. Computed mode 11 strain energy release rate 
distribution across the width of a SLB specimen for different models. 
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Figure 51. Computed moe III strain energy release rate 
distribution across the width of a SLB specimen for different models. 
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